#pragma once

#include "loirey_BasicStructure.h"
#include "loirey_MyMath.h"
#include "loirey_VectorMath.h"

namespace loirey
{

	namespace fft
	{

		extern int mUnitRoots;
		extern int nUnitRoots;
		extern CSimpleTypeArray<CComplexNumber> vUnitRoots;

		void ensure_unit_roots(int M);
		void set_complex_unit(int M, int k, CComplexNumber& Target);

		bool fft_Small_M_Only(
			int M, bool fInverse,
			const CComplexNumber* qSrc,
			CComplexNumber* pDst);

		void fft_Radix2(
			int M, bool fInverse,
			const CComplexNumber* qSrc,
			CComplexNumber* pDst);

		void fft_SplitRadix(
			int M, bool fInverse,
			const CComplexNumber* qSrc,
			CComplexNumber* pDst);

		void basic_perform_fft(
			int M, bool fInverse,
			const CComplexNumber* qSrc,
			CComplexNumber* pDst);

		const int BORDER_MODE_Border = 0;
		const int BORDER_MODE_Periodic = 1;
		const int BORDER_MODE_Zero = 2;

		const double COST_COEFFICIENT_FFT = 2.0;
		const double COST_COEFFICIENT_Convolution = COST_COEFFICIENT_FFT * 3.0;
		const double COST_COEFFICIENT_MultiplyPolynomial = COST_COEFFICIENT_FFT * 3.0;

		template<typename T>
		bool basic_format_data_1D(
			const T* qSrcData, int Size, int Step,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CComplexNumber* pDstData, int M, int N)
		{
			int i, tp;
			T td;

			if (BorderMode == BORDER_MODE_Border)
			{
				for ( i = 0; i < N; ++i ) pDstData[i].rel = pDstData[i].img = 0.0;
				for ( i = HeadIndex; i < TailIndex; ++i )
				{
					td = qSrcData[(i<0 ? 0 : (i>=Size ? Size-1 : i)) * Step];
					tp = ((fInverseData ? CenterIndex-i : i-CenterIndex) % N + N) % N;
					pDstData[tp] += td;
				}
			}
			else if (BorderMode == BORDER_MODE_Periodic)
			{
				for ( i = 0; i < N; ++i ) pDstData[i].rel = pDstData[i].img = 0.0;
				for ( i = HeadIndex; i < TailIndex; ++i )
				{
					td = qSrcData[(i % Size + Size) % Size * Step];
					tp = ((fInverseData ? CenterIndex-i : i-CenterIndex) % N + N) % N;
					pDstData[tp] += td;
				}
			}
			else if (BorderMode == BORDER_MODE_Zero)
			{
				for ( i = 0; i < N; ++i ) pDstData[i].rel = pDstData[i].img = 0.0;
				HeadIndex = (HeadIndex<0 ? 0 : HeadIndex);
				TailIndex = (TailIndex>Size ? Size : TailIndex);
				for ( i = HeadIndex; i < TailIndex; ++i )
				{
					td = ((i<0 || i>=Size) ? ((T)0.0) : qSrcData[i * Step]);
					tp = ((fInverseData ? CenterIndex-i : i-CenterIndex) % N + N) % N;
					pDstData[tp] += td;
				}
			}
			else
				return false;

			return true;
		}

		template<typename T>
		int basic_format_data_1D(
			const T* qSrcData, int Size, int Step,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CSimpleTypeArray<CComplexNumber>& vDstData, int MinLength)
		{
			int ret, N;
			for ( ret=0,N=1; N<MinLength; ++ret ) N <<= 1;
			vDstData.Resize(N);

			bool ok = basic_format_data_1D(
				qSrcData, Size, Step,
				HeadIndex, TailIndex, BorderMode,
				CenterIndex, fInverseData,
				vDstData.pElement, ret, N);
			if (!ok) ret = -1;

			return ret;
		}

		template<typename T>
		bool Perform_FFT_1D(
			const T* qData, int Size, int Step,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CComplexNumber* pResult, int M, int N,
			bool fInverseFFT)
		{
			CSimpleTypeArray<CComplexNumber> vData;
			vData.Resize(N);
			bool ret = basic_format_data_1D(
				qData, Size, Step,
				HeadIndex, TailIndex, BorderMode,
				CenterIndex, fInverseData,
				vData.pElement, M, N);
			if (ret)
				basic_perform_fft(M, fInverseFFT, vData.pElement, pResult);
			return ret;
		}

		bool Perform_FFT_1D(
			const CComplexNumber* qData, int Size, int Step,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CComplexNumber* pResult, int M, int N,
			bool fInverseFFT);

		template<typename T>
		int Perform_FFT_1D(
			const T* qData, int Size, int Step,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CSimpleTypeArray<CComplexNumber>& vResult, int MinLength,
			bool fInverseFFT)
		{
			int ret;
			int N;

			for ( ret=0,N=1; N<MinLength; ++ret ) N <<= 1;
			vResult.Resize(N);

			bool ok = Perform_FFT_1D(
				qData, Size, Step,
				HeadIndex, TailIndex, BorderMode,
				CenterIndex, fInverseData,
				vResult.pElement, ret, N,
				fInverseFFT);
			if (!ok) ret = -1;

			return ret;
		}

		template<typename T>
		int Perform_FFT_1D(
			const CSimpleTypeArray<T>& vData,
			int HeadIndex, int TailIndex, int BorderMode,
			int CenterIndex, bool fInverseData,
			CSimpleTypeArray<CComplexNumber>& vResult, int MinLength,
			bool fInverseFFT)
		{
			Perform_FFT_1D(
				vData.pElement, (int)vData.Size, 1,
				HeadIndex, TailIndex, BorderMode,
				CenterIndex, fInverseData,
				vResult, MinLength,
				fInverseFFT);
		}

		template<typename T>
		int Perform_FFT_1D(
			const CSimpleTypeArray<T>& vData,
			int CenterIndex, bool fInverseData,
			CSimpleTypeArray<CComplexNumber>& vResult, int MinLength,
			bool fInverseFFT)
		{
			Perform_FFT_1D(
				vData.pElement, (int)vData.Size, 1,
				0, (int)vData.Size, BORDER_MODE_Zero,
				CenterIndex, fInverseData,
				vResult, MinLength,
				fInverseFFT);
		}

		template<typename T>
		void Perform_FFT_2D(
			const T* p_matData, int W, int H,
			int xHead, int yHead, int xTail, int yTail, int BorderMode,
			int xCenter, int yCenter, bool fInverseData,
			CComplexNumber* pResult,
			int MW, int MH, int NW, int NH,
			bool fInverseFFT)
		{
			int ix, iy;

			CSimpleTypeArray<CComplexNumber> matCol;
			matCol.Resize(NH * W);

			for ( ix = 0; ix < W; ++ix )
			{
				Perform_FFT_1D(
					p_matData+ix, H, W,
					yHead, yTail, BorderMode,
					yCenter, fInverseData,
					matCol.pElement+NH*ix, MH, NH,
					fInverseFFT);
			}

			for ( iy = 0; iy < NH; ++iy )
			{
				Perform_FFT_1D(
					matCol.pElement+iy, W, NH,
					xHead, xTail, BorderMode,
					xCenter, fInverseData,
					pResult+NW*iy, MW, NW,
					fInverseFFT);
			}
		}

		template<typename T>
		void Perform_FFT_2D(
			const T* p_matData, int W, int H,
			int xHead, int yHead, int xTail, int yTail, int BorderMode,
			int xCenter, int yCenter, bool fInverseData,
			CSimpleTypeArray<CComplexNumber>& vResult, int MinW, int MinH,
			int& MW, int& MH, int& NW, int& NH,
			bool fInverseFFT)
		{
			for ( MW=0,NW=1; NW<MinW; ++MW ) NW <<= 1;
			for ( MH=0,NH=1; NH<MinH; ++MH ) NH <<= 1;
			vResult.Resize(NW * NH);

			Perform_FFT_2D(
				p_matData, W, H,
				xHead, yHead, xTail, yTail, BorderMode,
				xCenter, yCenter, fInverseData,
				vResult.pElement,
				MW, MH, NW, NH,
				fInverseFFT);
		}

		double CostConvolution_1D(
			int ResultHeadIndex, int ResultTailIndex,
			bool fExpandDataHead, bool fExpandDataTail,
			int MaskLength, int MaskCenterIndex);

		template<typename TData, typename TMask>
		void Convolution_1D(
			const TData* qData, int DataLength, int DataStep,
			int BorderMode,
			int ResultHeadIndex, int ResultTailIndex,
			bool fExpandDataHead, bool fExpandDataTail,
			const TMask* qMask, int MaskLength, int MaskStep,
			int MaskCenterIndex,
			CSimpleTypeArray<CComplexNumber>& vResult)
		{
			CSimpleTypeArray<CComplexNumber> vData;
			CSimpleTypeArray<CComplexNumber> vMask;
			int M, N, i;
			int DataHeadIndex = ResultHeadIndex - (fExpandDataHead ? MaskCenterIndex : 0);
			int DataTailIndex = ResultTailIndex + (fExpandDataTail ? MaskLength-1-MaskCenterIndex : 0);
			i = DataTailIndex - DataHeadIndex;
			for ( M=0,N=1; N<i; ++M ) N <<= 1;
			vData.Resize(N); vMask.Resize(N); vResult.Resize(N);

			Perform_FFT_1D(
				qData, DataLength, DataStep,
				DataHeadIndex, DataTailIndex, BorderMode,
				ResultHeadIndex, false,
				vData.pElement, M, N,
				false);
			Perform_FFT_1D(
				qMask, MaskLength, MaskStep,
				0, MaskLength, BORDER_MODE_Zero,
				MaskCenterIndex, true,
				vMask.pElement, M, N,
				false);
			for ( i = 0; i < N; ++i )
				vData.pElement[i] *= vMask.pElement[i];
			Perform_FFT_1D(
				vData.pElement, N, 1,
				0, N, BORDER_MODE_Zero,
				0, false,
				vResult.pElement, M, N,
				true);
			vResult.Resize(ResultTailIndex - ResultHeadIndex);
		}

		template<typename TData, typename TMask>
		void Convolution_1D(
			const TData* qData, int DataLength, int DataStep,
			int BorderMode,
			int ResultHeadIndex, int ResultTailIndex,
			bool fExpandDataHead, bool fExpandDataTail,
			const TMask* qMask, int MaskLength, int MaskStep,
			int MaskCenterIndex,
			CComplexNumber* pResult)
		{
			if (ResultHeadIndex >= ResultTailIndex) return;

			CSimpleTypeArray<CComplexNumber> vTemp;
			Convolution_1D(
				qData, DataLength, DataStep,
				BorderMode,
				ResultHeadIndex, ResultTailIndex,
				fExpandDataHead, fExpandDataTail,
				qMask, MaskLength, MaskStep,
				MaskCenterIndex,
				vTemp);

			DataLength = ResultTailIndex - ResultHeadIndex;
			int i;
			for ( i = 0; i < DataLength; ++i )
				pResult[i] = vTemp.pElement[i];
		}

		template<typename TData, typename TMask, typename TResult>
		void Convolution_1D(
			const TData* qData, int DataLength, int DataStep,
			int BorderMode,
			int ResultHeadIndex, int ResultTailIndex,
			bool fExpandDataHead, bool fExpandDataTail,
			const TMask* qMask, int MaskLength, int MaskStep,
			int MaskCenterIndex,
			TResult* pResult)
		{
			if (ResultHeadIndex >= ResultTailIndex) return;

			CSimpleTypeArray<CComplexNumber> vTemp;
			Convolution_1D(
				qData, DataLength, DataStep,
				BorderMode,
				ResultHeadIndex, ResultTailIndex,
				fExpandDataHead, fExpandDataTail,
				qMask, MaskLength, MaskStep,
				MaskCenterIndex,
				vTemp);

			DataLength = ResultTailIndex - ResultHeadIndex;
			int i;
			for ( i = 0; i < DataLength; ++i )
				pResult[i] = CMyMath::TranslateNumber_Round(vTemp.pElement[i].rel, (TResult)0);
		}

		template<typename TData, typename TMask, typename TResult>
		void Convolution_1D(
			const TData* qData, int DataLength, int DataStep,
			int BorderMode,
			int ResultHeadIndex, int ResultTailIndex,
			bool fExpandDataHead, bool fExpandDataTail,
			const TMask* qMask, int MaskLength, int MaskStep,
			int MaskCenterIndex,
			CSimpleTypeArray<TResult>& vResult)
		{
			if (ResultHeadIndex >= ResultTailIndex) return;

			vResult.Resize(ResultTailIndex - ResultHeadIndex);
			Convolution_1D(
				qData, DataLength, DataStep,
				BorderMode,
				ResultHeadIndex, ResultTailIndex,
				fExpandDataHead, fExpandDataTail,
				qMask, MaskLength, MaskStep,
				MaskCenterIndex,
				vResult.pElement);
		}

		double CostMultiplyPolynomial(
			int LengthA, int LengthB);

		template<typename TA, typename TB>
		void MultiplyPolynomial(
			const TA* qA, int LengthA, int StepA,
			const TB* qB, int LengthB, int StepB,
			CSimpleTypeArray<CComplexNumber>& vResult)
		{
			CSimpleTypeArray<CComplexNumber> vA;
			CSimpleTypeArray<CComplexNumber> vB;
			int M, N, i;
			i = LengthA + LengthB - 1;
			for ( M=0,N=1; N<i; ++M ) N <<= 1;
			vA.Resize(N); vB.Resize(N); vResult.Resize(N);

			Perform_FFT_1D(
				qA, LengthA, StepA,
				0, LengthA, BORDER_MODE_Zero,
				0, false,
				vA.pElement, M, N,
				false);
			Perform_FFT_1D(
				qB, LengthB, StepB,
				0, LengthB, BORDER_MODE_Zero,
				0, false,
				vB.pElement, M, N,
				false);
			for ( i = 0; i < N; ++i )
				vA.pElement[i] *= vB.pElement[i];
			Perform_FFT_1D(
				vA.pElement, N, 1,
				0, N, BORDER_MODE_Zero,
				0, false,
				vResult.pElement, M, N,
				true);
			vResult.Resize(LengthA+LengthB-1);
		}

		template<typename TA, typename TB>
		void MultiplyPolynomial(
			const TA* qA, int LengthA, int StepA,
			const TB* qB, int LengthB, int StepB,
			CComplexNumber* pResult)
		{
			CSimpleTypeArray<CComplexNumber> vTemp;
			MultiplyPolynomial(
				qA, LengthA, StepA,
				qB, LengthB, StepB,
				vTemp);

			int i;
			for ( i = 0; i < (int)vTemp.Size; ++i )
				pResult[i] = vTemp.pElement[i];
		}

		template<typename TA, typename TB, typename TResult>
		void MultiplyPolynomial(
			const TA* qA, int LengthA, int StepA,
			const TB* qB, int LengthB, int StepB,
			TResult* pResult)
		{
			CSimpleTypeArray<CComplexNumber> vTemp;
			MultiplyPolynomial(
				qA, LengthA, StepA,
				qB, LengthB, StepB,
				vTemp);

			int i;
			for ( i = 0; i < (int)vTemp.Size; ++i )
				pResult[i] = CMyMath::TranslateNumber_Round(vTemp.pElement[i].rel, (TResult)0);
		}

		template<typename TA, typename TB, typename TResult>
		void MultiplyPolynomial(
			const TA* qA, int LengthA, int StepA,
			const TB* qB, int LengthB, int StepB,
			CSimpleTypeArray<TResult>& vResult)
		{
			vResult.Resize(LengthA+LengthB-1);
			MultiplyPolynomial(
				qA, LengthA, StepA,
				qB, LengthB, StepB,
				vResult.pElement);
		}

		double CostConvolution_2D(
			int Result_xHead, int Result_yHead, int Result_xTail, int Result_yTail,
			bool fExpand_xHead, bool fExpand_yHead, bool fExpand_xTail, bool fExpand_yTail,
			int MaskW, int MaskH,
			int Mask_xCenter, int Mask_yCenter);

		template<typename TData, typename TMask>
		void Convolution_2D(
			const TData* p_matData, int DataW, int DataH,
			int BorderMode,
			int Result_xHead, int Result_yHead, int Result_xTail, int Result_yTail,
			bool fExpand_xHead, bool fExpand_yHead, bool fExpand_xTail, bool fExpand_yTail,
			const TMask* p_matMask, int MaskW, int MaskH,
			int Mask_xCenter, int Mask_yCenter,
			CSimpleTypeArray<CComplexNumber>& matResult)
		{
			CSimpleTypeArray<CComplexNumber> matData;
			CSimpleTypeArray<CComplexNumber> matMask;
			int MW, MH, NW, NH, i;
			int RW = Result_xTail - Result_xHead;
			int RH = Result_yTail - Result_yHead;
			int Data_xHead = Result_xHead - (fExpand_xHead ? Mask_xCenter : 0);
			int Data_yHead = Result_yHead - (fExpand_yHead ? Mask_yCenter : 0);
			int Data_xTail = Result_xTail + (fExpand_xTail ? MaskW-1-Mask_xCenter : 0);
			int Data_yTail = Result_yTail + (fExpand_yTail ? MaskH-1-Mask_yCenter : 0);
			for ( MW=0,NW=1; NW<Data_xTail-Data_xHead; ++MW ) NW <<= 1;
			for ( MH=0,NH=1; NH<Data_yTail-Data_yHead; ++MH ) NH <<= 1;
			matData.Resize(NW*NH); matMask.Resize(NW*NH); matResult.Resize(RW*RH);

			Perform_FFT_2D(
				p_matData, DataW, DataH,
				Data_xHead, Data_yHead, Data_xTail, Data_yTail, BorderMode,
				Result_xHead, Result_yHead, false,
				matData.pElement, MW, MH, NW, NH,
				false);
			Perform_FFT_2D(
				p_matMask, MaskW, MaskH,
				0, 0, MaskW, MaskH, BORDER_MODE_Zero,
				Mask_xCenter, Mask_yCenter, true,
				matMask.pElement, MW, MH, NW, NH,
				false);
			for ( i = 0; i < (int)matData.Size; ++i )
				matData.pElement[i] *= matMask.pElement[i];
			Perform_FFT_2D(
				matData.pElement, NW, NH,
				0, 0, NW, NH, BORDER_MODE_Zero,
				0, 0, false,
				matMask.pElement, MW, MH, NW, NH,
				true);
			for ( i = 0; i < (int)matResult.Size; ++i )
				matResult.pElement[i] = matMask.pElement[i%RW+i/RW*NW];
		}

		template<typename TData, typename TMask>
		void Convolution_2D(
			const TData* p_matData, int DataW, int DataH,
			int BorderMode,
			int Result_xHead, int Result_yHead, int Result_xTail, int Result_yTail,
			bool fExpand_xHead, bool fExpand_yHead, bool fExpand_xTail, bool fExpand_yTail,
			const TMask* p_matMask, int MaskW, int MaskH,
			int Mask_xCenter, int Mask_yCenter,
			CComplexNumber* p_matResult)
		{
			CSimpleTypeArray<CComplexNumber> matTemp;
			Convolution_2D(
				p_matData, DataW, DataH,
				BorderMode,
				Result_xHead, Result_yHead, Result_xTail, Result_yTail,
				fExpand_xHead, fExpand_yHead, fExpand_xTail, fExpand_yTail,
				p_matMask, MaskW, MaskH,
				Mask_xCenter, Mask_yCenter,
				matTemp);

			int i;
			for ( i = 0; i < (int)matTemp.Size; ++i )
				p_matResult[i] = matTemp.pElement[i];
		}

		template<typename TData, typename TMask, typename TResult>
		void Convolution_2D(
			const TData* p_matData, int DataW, int DataH,
			int BorderMode,
			int Result_xHead, int Result_yHead, int Result_xTail, int Result_yTail,
			bool fExpand_xHead, bool fExpand_yHead, bool fExpand_xTail, bool fExpand_yTail,
			const TMask* p_matMask, int MaskW, int MaskH,
			int Mask_xCenter, int Mask_yCenter,
			TResult* p_matResult)
		{
			CSimpleTypeArray<CComplexNumber> matTemp;
			Convolution_2D(
				p_matData, DataW, DataH,
				BorderMode,
				Result_xHead, Result_yHead, Result_xTail, Result_yTail,
				fExpand_xHead, fExpand_yHead, fExpand_xTail, fExpand_yTail,
				p_matMask, MaskW, MaskH,
				Mask_xCenter, Mask_yCenter,
				matTemp);

			int i;
			for ( i = 0; i < (int)matTemp.Size; ++i )
				p_matResult[i] = CMyMath::TranslateNumber_Round(matTemp.pElement[i].rel, (TResult)0);
		}

		template<typename TData, typename TMask, typename TResult>
		void Convolution_2D(
			const TData* p_matData, int DataW, int DataH,
			int BorderMode,
			int Result_xHead, int Result_yHead, int Result_xTail, int Result_yTail,
			bool fExpand_xHead, bool fExpand_yHead, bool fExpand_xTail, bool fExpand_yTail,
			const TMask* p_matMask, int MaskW, int MaskH,
			int Mask_xCenter, int Mask_yCenter,
			CSimpleTypeArray<TResult>& matResult)
		{
			int RW = Result_xTail - Result_xHead;
			int RH = Result_yTail - Result_yHead;
			matResult.Resize(RW*RH);
			Convolution_2D(
				p_matData, DataW, DataH,
				BorderMode,
				Result_xHead, Result_yHead, Result_xTail, Result_yTail,
				fExpand_xHead, fExpand_yHead, fExpand_xTail, fExpand_yTail,
				p_matMask, MaskW, MaskH,
				Mask_xCenter, Mask_yCenter,
				matResult.pElement);
		}

	}

}

